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This is the first of a series of papers describing a numerical implementation 
of the conformally rescaled Einstein equation, an implementation designed 
to calculate asymptotically flat spacetimes, especially spacetimes containing 
black holes. 

In the present paper we derive the new first order time evolution equations to 
be used in the scheme. These time evolution equations can either be written 
in symmetric hyperbolic or in flux-conservative form. Since the conformally 
rescaled Einstein equation, also called the conformal field equations, formally 
allow us to place the grid boundaries outside the physical spacetime, we can 
modify the equations near the grid boundaries and get a consistent and stable 
discretisation. Even if we calculate spacetimes containing black holes, there is 
no need for introducing artifical boundaries in the physical spacetime, which 
then would complicate, influence, or even exclude the computation of certain 
spacetime regions. 

I. INTRODUCTION 



Since the capacity of computers has been rapidly growing over the last years, numerical 
experiments on spacetimes without symmetries are now in the reach of the available hardware 
resources. Numerical relativity is on the way to take the step which has already been taken 
in other fields of physics: Numerical experiments supplement, or sometimes even replace, 
real experiments. 

When designing the computer programs for the numerical experiments it is highly important 
to be aware of the following significant difference to other fields of physics: One of the 
major methods in developing codes for numerical experiments has been the comparison and 
testing with real experiments. In numerical relativity however, such possibilities are very 
rare and limited; appropriate experiments are not available, comparison with observations 
of astrophysical events requires the knowledge of many unknown parameters, and the known 
exact solutions presumably do not represent the full spectrum of properties of solutions of 
the non-linear Einstein equation. 

The author, therefore, believes, that reliability of the results of numerical experiments in 



relativity can only be ensured by mathematical rigor of the implemented algorithms, i. e. 
their vanishing grid size limit should be consistent with the Einstein equation. 

In this series of papers we describe a numerical implementation of the conformal field 
equations [|l],@. In this implementation we do not use any approximation besides the ideal- 
isation of isolated systems and the replacement of the partial differential equations by their 
discretised versions. The approach described below is designed to be able to calculate global 
as well as local properties of spacetimes arising from data of arbitrary strength. Important 
issues are e. g. testing the validity range of approximations, calculating the long-time be- 
haviour of fields, and investigating the nature of singularities. For spherically symmetric 
models with scalar fields the attainability of these goals has been demonstrated in 
The approach described is so powerful that we can use the same time evolution scheme to 
calculate scenarios like the decay of weak gravitational waves to Minkowski space or the 
merging of black holes. 

The use of the conformal field equations allows us to extend the calculation to null infinity. 
Moreover, we can formally even extend the spacetime through null infinity and place the grid 
boundaries in a region which is causally disconnected from the physical spacetime. Then 
boundary effects, like a violation of the constraints or creation of spurious backscattering, 
cannot propagate into the physical spacetime. As null infinity is part of the grid, we can 
calculate anything which is suggested by a Penrose diagram, including quantities which are 
only defined in the limit of going to infinity, like gravitational radiation. 

Since an exhaustive description of the new ingredients needed for this (n+l)-dimensional 
(n<3) implementation of the conformal approach is too long for a single paper and requires 
techniques from mathematical relativity as well as numerical mathematics, we have split it 
into a series of papers. In this first part of the series we describe the foundation of the time 
evolution part of the code, namely the geometrical background, the equations used, and the 
physical scenarios which will be treated in detail in the following articles. In particular, we 
develop the mathematical basis and give a preview over the result which can be expected 
from our numerical approach. 

Details about the initial data problem, the discretisation used, the tests, and more will be 
found in the following parts of the series. 

In section |I| the conformal field equations are written in first order form and split into 
symmetric hyperbolic time evolution equations and constraints. When deriving this new 
form, which could also be written in flux conservative form, we proceed in analogy to ||, 
a paper dealing with the Einstein equations for vacuum in physical spacetime. Also, these 
properties of the conformal field equations are recalled which are essential for the follow- 
ing. Having the conventions available we discuss the hyperboloidal initial value problem 



in section |T|. There we also explain the special role of the conformal factor, special not 
as a variable of the system, but special when the physical spacetime is reconstructed from 
the variables of the conformal field equation. Section |IV| describes possible grid boundary 
treatments and proves their numerical stability. 

In the last section we outline the power of the conformal approach by describing why we can 
in principle calculate the complete future evolution of the initial data, including null infin- 
ity. Since conformal rescaling preserves the null cone structure of a spacetime, the conformal 
approach is especially tuned for analysing causal structure and gravitational radiation. 
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II. A SYMMETRIC HYPERBOLIC FIRST ORDER SYSTEM FOR THE 
CONFORMAL FIELD EQUATIONS 

We start the derivation of the first order equations by giving the conformal field equations 
in our notation and recalling some of their properties, whose proofs can be found in HQ. 



A. The conformal field equations and the reconstruction of physical spacetime 

In the conformal approach we solve the initial value problem for the conformal field 
equations, which read 

V a R bc - V b R ac + — {{V a R) g bc - {V b R) g ac ) + 2 (V d fi) d abc d = 0, (la) 

Vdd abc d = 0, (lb) 
V a V b n a + ^R ab n-^V a V a Qg ab = 0, (lc) 

\ V Q (v b V b fi) + l - R ab V b Q + ^Q V a R +^V a QR = 0, (Id) 

R 

fld ab d + (g ca R b d - g cb R a d - g d a R bc + g d b R ac )/2 + (g ca g b d - g cb g d )- = R abc d , (le) 
and 

tt 2 R + 6QV a V a n- 12(V a fi) (V fl fi) = 0. (If) 

Here the tensor g ab is a Lorentzian metric on a manifold M, R ab a symmetric-traceless tensor 
field, d abc d a tensor field with the symmetries of the Weyl tensor, fl a scalar, R abc d the differ- 
ential expression of the Riemann tensor in terms of the metric g ab , R a prescribed function 
on M, and V a denotes the derivative operator associated to g ab . 

For any solution (g ab , R ab , d abc d , fl) of system ([!]) the function R turns out to be the Ricci 

scalar, R ab the traceless part of the Ricci tensor, and fld abc d the Weyl tensor of g ab . 

Any solution of (p!a|-[ldD solves (flfj) on a connected domain if (0) holds at one point in the 



domain. It is, therefore, sufficient to enforce ( |TI]) as a constraint on the initial data. In the 
following this requirement will no longer be explicitly mentioned. 

We will arrange our calculation such that (M, g ab ) will be a globally hyperbolic manifold 
sliced by the spacelike hypersurfaces E t . By solving ([]]) on M we obtain a solution of the 
Einstein equation on the physical spacetime M := {x G M | Q(x) > 0}, since for any suf- 
ficiently smooth solution of ([Ia|-[lcJ) the metric g ab := Q~ 2 g ab is a solution of the Einstein 



equation on M. 

It is sufficient for our purposes to only consider setups for which the physical part 
E to := {x G S to | ft(x) > 0} of the initial hypersurface T, to is relatively compact. We 
define S t := {x G T, t \ Q(x) = 0}. 

The part of the boundary of M with V a fi \n=o^ is denoted J . All geodesies of g ab which 
end at J are null geodesies. Since they intersect J at infinite affine parameter with respect 
to g ab , J is infinitely far away with respect to g ab . J represents, therefore, null infinity, 
and quantities like gravitational radiation, which is extracted most effectively by taking the 
limit of going to infinity, can be determined by reading off the values of certain variables at 
J. 
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As can be seen immediately from equation flip, J is a null hypersurface. When solving the 
initial value problem for data given on the spacelike hypersurface S to , we will, unless men- 
tioned otherwise, require V a f2 to be future directed everywhere on S := {P G | Q(P) = 
0}. Or stated more interpretatively, we choose the initial hypersurface in such a way that 
it intersects future null infinity and not past null infinity. Then £ <0 : = M D S to is a Cauchy 
surface for the part of M which lies in the future of S io , i. e. M n P(E to ) = X>(M n £t ), 
where P denotes the future domain of dependence.^ This property guarantees that at the 
points of {a; | Q(x) > 0} the solution of the initial value problem is independent of the values 
of the variables at the points in {x\ Q(x) < 0}. Hence, by placing the grid boundaries outside 
physical spacetime, namely into the set {a;| Q(x) < 0}, we gain a decoupling of the boundary 
treatment from the calculation of physical spacetime. 

If J = R x S consists of (disconnected) pieces with topology Rx S 2 , the spacetime is called 
(weakly) asymptotically flat ||. The gravitational wave scenarios usually treated fall into 
this class. The conformal field equations also allow us to study data with other topologies 
of S, and in fact many aspects of gravitational radiation can also be studied on such space- 
times. We will see in subsection [111 B 2[ that in numerical studies spacetimes with toroidal 
sections of J may provide advantages over those with J' = R x S 2 , especially in the cases 
where we reduce the space dimensions by assuming symmetries. To have a term available 
which refers to both cases we call spacetimes asymptotically regular ,0 if they arise from 
regular data for the conformal field equations with V a f2 future directed everywhere on S to . 
The conformal and the physical metric are related by a rescaling which is essentially arbri- 
trary, as two solutions (M, gab, ft) and (M, g a b,&) with (g a b, fl) = (9 2 g a b,9Q) and a positive 
function 9 describe the same physical spacetime. Under the rescaling 9 the Ricci scalar 
R changes. To obtain equations with nice properties we do not prescribe Q but R. The 
prescribed function R is therefore a gauge source function for the rescaling gauge. 
The equations (|TaD , (|lb| ) and (|ld| ) are third order equations for the metric g a b and the con- 
formal factor Q. To be able to apply standard theorems which imply well-posedness of the 
initial value problem, we introduce appropriate gauge conditions and formulate the equa- 
tions (|l|) as a symmetric hyperbolic first order system. There are various possibilities to do 
this. In all variants gauge source functions must be specified to make the solution of the 
initial value problem unique. In the variants which use the spinor ||[7j or the frame formal- 
ism [p|,|T0f the coordinate gauge freedom is fixed by 10 gauge source functions (4 degrees of 
freedom for the coordinates, 6 for the tetrad). In the variant to be described here, which 
uses a 3-tensor formalism, the coordinate gauge freedom is fixed by the three components 
of the shift vector and a function which is closely related to the lapse function. 
The possible implementations of the conformal field equations differ slightly in the number 
of variables, the gauge conditions, and the type of non-linearity. Our propagation equations 
allow us to write them in flux conservativ form. Whether one of the implementations has 
significant advantages over the others in numerical calculations remains to be seen. A nu- 



1 V(S) := {P € M | Every past inextendible causal curve through P intersects S} 

2 regular with respect to the existence of null infinity, the spacetime may nevertheless contain 
singularities. 
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merical comparison with [|7|,[TTJ is planned (see also subsection |IIIB2| ). 
Many quantities of physical interest, like gravitational radiation or tidal forces (geodesic 
deviation), can be expressed in terms of the components of the conformal Weyl tensor d a b c d 
or the physical Riemann tensor. As the systems just mentioned contain the components of 
the conformal curvature as variables, those quantities can be calculated by pure algebra and 
their error can be estimated by a convergence analysis. 



B. The 3+1 split of the conformal field equations 

1. Equations for the lower order quantities 

Proceeding in analogy to [[5] we introduce the extrinsic curvature k a b as new variable. It 
is related to the 3-metric0 hob and the hypersurface normal 3 n a by the equation 

C n h ab - 2k ab = 0. (2a) 

The connection T a b c relates to the curvature R a bc d by 

— dcX d bc + dbT d ac + T e ca T d be — T e c bT d ae — R a bc = 0. (2b) 

Equation (|lc| ) and (|Td"D, the equations for the conformal factor Q, are written in first order 
form as follows, 

V a fi - (4 U a = 0, (2c) 
V a (4 to 6 + -iL,O-a;<? a6 = 0, (2d) 

with cu = iV a V a fi, and 

V a u + ~ R ah {i Xl b + 1 Q V a R + ^ (4) fi a R = 0. (2e) 

In the equation for the tracefree Ricci tensor R a b we substitute the derivatives of fl, 

V [a R b]c + ^(y la R) g b]c + i4 ^ d d abc d = 0. (2f) 

The equation for the conformal Weyl tensor d a b c d , 

V d d abc d = 0, (2g) 

is kept. 

We now make a 3+1 split of the system (2) by a straightforward extension of the treatment 
given in 0. For completeness the essential definitions are repeated but for further clarifica- 
tion the reader is advised to consult this reference. 

The globally hyperbolic manifold (M,g a b) is sliced with spacelike hypersurfaces S t . The 



3 definition follows 
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parameter t then provides a natural time coordinate and we can introduce a timelike vector 
field t a related to t by t a V a t = 1. Denoting the unit normal of the hypersurface E t by n a , 
n a and t a are related by 

n a = ^(t a -N a ), (3) 

with A" the lapse function and iV a the shift vector. The metric g ab induces a 3-metric h ab 
on each S t by 

gab = h ab - n a n b . (4) 

We define 

a a := n b V b n a . (5) 

The object a a = ha ® bN is purely spatial. 4 Later the coordinate gauge freedom is fixed by 
specifying the three components of the shift vector N a and the function q := ln(^=) as free 
functions of space and time. The following holds: 

a a = h a b d b q + 7\ fe . (6) 

7 a b c is the 3-connection for the 3-metric h ab and is used as a variable of the first order system. 
The 4-connection T a bc is decomposed into 3-parts by means of the hypersurface normal, 

T\ c = {1,1 ' iy r a bc - n a ^T bc - (1 ' 0,1) r a c n b - ^T a b n c 

+ (1 '°'°T a n b n c + n a (0 '°' 1 T b n c + n a (0 '°' 1 T C n b - n a (0 '°'°T n b n c , (7) 

where a (-'°'-) in front of a quantity indicates contraction with n a and a means 
projected with h a b , e.g. (1,0 ' 1) T o (; := h a a-n b h c e T d be . A calculation shows 

(8a) 

(8b) 

(8c) 
(8d) 

(8e) 
(8f) 

For the decomposition of the curvature variables we get 

Rab = (1A) R ab - n a ^R b - ™R*n b + n a n b (0 ' 0) R (9) 

and 



(0,0, 0)p 






AT 
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4 can be calculated on Y± t by only knowing the values of N, N a , and h a b on E< 
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dabcd — IdbEac — IdaEbc — ld>E a b + l Ca E bd 

- n b B ae ^e e dc + n a B b} <% e dc - n d B cf "\ e ba + n c B df < :i k\ a , (10) 
with l ab = Kb + n a n b , {3) e bcd = n a e abcd , E a a = 0, and B a a = 0. As R ab is tracefree, 

For the derivatives of f2 we substitute 

n ■= n a V a n = n a{i \l a (11) 

and 

n a := h a b v b n = h a b ^ b . (12) 

The set of 3-variables is now complete, it is (h a b, k a b, J a bc, (0,1> Ra, (1 ' 1} Rab, E ab , B ab , Q, fi , 
Q a , u). 

2. The split into time evolution equations and constraints 

Contracting all indices of the equations (2) with n a respectively h b and using the new 
variables we obtain time evolution equations (equations with d t ) and constraints (equations 
without d t ). Although this calculation is lengthy, it is straightforward. To perform it 
the author has used a MATHTENSOR program for doing 3+1 decompositions written by 
himself. Then, after removing the gauge freedom by fixing gauge source functions and adding 
constraints to the time evolution equations to achieve symmetric hyperbolicity, we obtain a 
complete set of symmetric hyperbolic time evolution equations. We write the equations in 
terms of the null quantities Af, the corresponding equations are obtained by setting the null 
quantities to 0: 

A" h ab = - £ n hab + 2k ab (13a) 
•Mc ab = - Cnkab + {3) ^ cl° ab + Ybc^ad + a a a b + k c c k ab - Y ' ab a c 

+ h a c h b d d d d c q - ^h ab - ^R c c h ab - 2QE ab (13b) 
A/" 7 a bc = - h ad Cra d bc + ^ahc ~ a-ahc + a c k ab + a b k ac + h da h b e hj : —d f d e N d 

- \h ac ^R b - hi ab ^R c + h bc ^R a + {: % c d nB db + (3) e ab d nB dc (13c) 
A/e ab = — C n E ab + \ (:i) e a cd ^ d B cb + \ {: % cd (:i V d B ca + a c < : % b d B da + a c < : % a d B db 

- h ab k cd E cd + \k h c E ca + h a c E cb - 2k c c E ab (13d) 
A/b ab — — CnB a b - \ (3) e a cd (3 V ' d E cb - \ ^e b cd (3 V d E ca + a c '\ b c d E da + a c { \ ac d E db 

- h ab k cd B cd + h b c B ca + b -k a c B cb - 2k c c B ab (13e) 

Mo,ift a = ~ Cn (ihl) Ra + ^b ^Ra" ~ \ ^ aR - k b b ^R a + a b ^R b + d a ^R b b (13f ) 

ab = - h bc C n ^R a c + (a V a ^R b - hi ab C n R + a b (0A) R a + a a ^R b 
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- k ab ^R c c - k cb ^R a c - 2 { \ cd b n c B da - 2Q E ab (13g) 
A/n = - + (13h) 

A/n = " AA) - w + a a Q a - 1 1 1 - 1 ^ (131) 

A/a a a = " £ A + a A + - ^ (M) ^a (13j) 

AC = - C n u - ^C n R - ^fio - y (M K + y CM) iC (13k) 

Note that for a 3-tensor t... a ..."' b "' the following equality holds: 

C.„l ...„... "'"" = — (Ct - C v ) t... a ...'" b '" + . . . + //'/...„..."""" a c + . . . 

The operator (3 V a is the covariant derivative preserving the 3-metric h ab . 
If one uses new variables for E ab and B ab as in [H] , this system is indeed symmetric hyperbolic. 
The equations ( |13a| - |13e] ) have the same principal part as the corresponding equations (6.4- 
6.8) in ||. The principal parts of ( |13a| ) and ( |13h| - |T3"k| ) are equal, the same is true for 
the principal parts of the subsystems (|13b , 13c|) and (|131| , |13gj ). Therefore the system (13 



possesses the same characteristics as the equations (6.4-6.8) in M, namely they are either 
null hypersurfaces with respect to g ab or they are timelike and tangent to n a or to the timelike 
cone {c ab t a t b = 0} where c ab = 4h ab — n a n b . Since all characteristics are within the light 
cone of g ab , the system is consistent with Einstein causality. 

The treatment of the coordinate gauge freedom (fixed by the functions q and N a ) is also 
discussed in more detail in . What is said there extends straightforwardly to the rescaling 
freedom fixed by the gauge source function R. 

For the numerical treatment the equation ( |13cj ) is multiplied by h ea and the h bc in ( |13g| ) 
is moved into the Lie derivative. Both manipulations are done to get times the identity 
matrix in front of the time derivative and the second is also done to obtain the 6 components 
of the symmetric tensor (1,1) -R a b instead of the 9 components of {1 ' 1} R a b as variables of the 
system. The two equations resulting from A/"(i,i& and A/"(i,i£ are identical only up to 
constraints and we take the arithmetic mean of the two null quantities. By this averaging 
some characteristics change but they remain within the light cone of g ab . 
The system (|13D can also be written in flux conservative form. The results we obtained 
by solving the flux conservative system with a rotated Richtmyer scheme were less accurate 
than those obtained by solving the quasilinear form (0) with the same scheme. Furthermore 
calculations on a simple non-linear model system showed that damping of grid modes is 
weaker for the flux-conservative form. Therefore, the quasilinear form of the equations 
seems to be better suited for our application than the flux conservative form and we refrain 
from presenting the flux conservative form of the system. 

The null quantities not appearing in the system ([H|) are the constraints: 

Nhabc = (3 *J a h bc (14a) 
A/k abc = - {3 ^ah c + (3 V b k ac + l -h ca ^R b - h cb ^R a - (3) e ab d QB dc (14b) 

AC abc ^ 'at be '^feT ac T ae~) be T be~) ac 
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k a k bc + k ac k b + Y2^ a ^ bc ^ ~^j l achb R 

- -h b d (M) iL + -h a d (M) £ 6c - -h ac ^R b d + -h a d ^R bc 
2 2 2 2 

- + M^a" + h d flE ac - h a d flE bc (14c) 
A/e a = — i3 V b E a b - ^e abc k bd B d c (14d) 
A/Ba = - (3 V fe S a b + (3) e abc ^£/ (14e) 

A/" ( o,i)r a6 = (3 V a ^R b - (3 V fe (cu) i? a + k b c ^R ca - k a c ^R cb + 2 { \ ab c fl d B d (14f ) 

- {0 ' 1] R b k ac + 2 (3) e ab d fi 5 dc - 2fl a E bc + 2ft b £ ac + 2h ca fl d E b d - 2h cb fl d E d (14g) 
A/" n «= - (3 V a fi + fi a (14h) 

Nn oa = - (3 V a fi + fcA - \fl ^Ra (141) 

A/h. at = - <3 VA + /l afe W + k ab fl - ^fi (la) i? a;) (I4j) 
K a = ~ (3 V a U - ^fl (3 V a R - ^fl a R + h ^Ra - l -tt h ^R ba (14k) 

The constraints are, just as the time evolution equations, regular on the whole manifold M 
including J . We abstain from proving the propagation of the constraints under the time 
evolution equations for the first order form of the conformal field equations derived. Since 
the propagation of the constraints is supported by the following circumstances, there is no 
serious doubt that they do propagate. First, the propagation has been proven for the forms 
of the conformal field equations given in fI^ , |T0| j7H . Second, for the first order form of the 



Einstein equation derived in ||, which is obtained from the form given by setting R = 0, 
(0,1) R a = 0, (1,1) i? a fe = 0, fl = 1 and hence d abc d = C abc d , the propagation of the constraints 
has been proven. And, third, under the numerical time evolution the constraints propagate 
modulo the order of convergence of the scheme used. 

To give initial data we have to find a solution of the constraints. With the exception of cases 
with high symmetry 0, it is not known yet how to solve the constraints (|l"4D directly for the 
variables / := (h ab , k ab , 7 a f, c , E ab , B ab , (0,1) i? a , {1,1) R ab , fl, fl , fl a , u). Nevertheless the existence 
of solutions of the constraints on S fl <9£ has been proven in [jnj for arbitrary topology of 
<9E. After solving an elliptic problem for (h ab , k ab , fl, flo), appropriate contractions of (|14"D 
can be used to calculate the remaining data. It is shown in JT3[, that, although divisions 



by fl are needed to calculate the curvature related quantities, these curvature quantities are 
smooth at <9S, where fl vanishes (= S), and provide, therefore, suitable data for the time 
evolution equations. Clearly, this division by fl is numerically a delicate issue. However, 



in [14| we will discuss a way of doing this which works for arbitrary S. 



III. THE ROLE OF THE CONFORMAL FACTOR 

In the preceeding section we have derived symmetric hyperbolic evolution equations of 
first order representation from the conformal field equations. Apart from the treatment of 
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the grid boundaries there are standard techniques available to obtain stable discretisations 
of symmetric hyperbolic systems. Due to properties distinguishing the conformal approach 
from approaches working in terms of physical spacetime the treatment of the outer grid 
boundary can be made trivial, as we will see in the next section. Therefore, it is interesting 
to understand what kind of physical scenarios can be dealt with by the conformal approach 
without complicating matters by introducing additional (inner) boundaries. 
We will see in this section that the class of scenarios meeting the requirements just mentioned 
is much larger than the class we would get under these requirements when working in terms 
of physical spacetime. The underlying reason is the role of the conformal factor Q. For the 
equations, and therefore for the discretisations, the conformal factor is just one of many 
variables which may assume any value. But the conformal factor is a very special variable 
when transforming back to physical spacetime by rescaling g ab \— > Q~ 2 g ab . Then the Q = 
level sets turn out to be boundaries of the physical spacetime. This double role of the 
conformal factor not only enlarges the class of physical problems which we can treat by 
standard techniques of numerical mathematics, it also determines the character of the initial 
value problem we solve. 

A. The hyperboloidal initial value problem 

Let f(t ,x) = (h ab ,k ab ,-f a bc ,E ab ,B ab , (0A) R a , {1A) R a c , fi, Oq, tt a , u) be a smooth solution 
of the constraints on an everywhere spacelike hypersurface E 4o , Q |s to > 0, with boundary 
<St , where Q vanishes. If V a fi |s is a future directed, non- vanishing null vector, we call 
(£t , f{to, %)) a hyperboloidal initial value problem. 

We can smoothly extend the hypersurface E to and the data f(t , x) beyond the boundary to 
obtain an extended hyperboloidal initial value problem (T, to , f(t , x)). Without loss 
of generality we assume that Q < on the extension. For an extended hyperboloidal initial 
value problem the physical part of the future of E to is independent of the formal extension 
of the initial data beyond S since M D 22(E) = T>{M n E). Therefore, we do not require the 
extended data to satisfy the constraints outside E to . 

B. Interpretation of various topologies of S 

In this subsection we classify various initial data configurations by finding spacetimes 
with the same topology of S. It may well be that the data in each class develop into 
spacetimes with very different structures in the large. To analyse this variety in detail is an 
interesting issue and it is a goal of the numerical experiments. Here, our purpose is to only 
describe a minimal collection of scenarios which we can treat with the conformal approach 
without introducing inner boundaries. 

1. Ss with spherical boundaries 

The simplest case of a compact three dimensional manifold S U E is shown in figure [I], 
namely the ball B 3 . 
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FIG. 1. Data evolving into a space- 
time with an asymptotic structure as 
Minkowski spacetime. The third di- 
mension is suppressed, S is a sphere. 




Grid boundary 



FIG. 2. Data evolving into a space- 
time with an asymptotic structure as 
Schwarzschild-Kruskal. The third di- 
mension is suppressed. The «S's are 
spheres, the corresponding event hori- 
zons are also shown. 



Its interior is diffeomorphic to R 3 . For a standard choice of the metric this is a hyper- 
boloidal slice in Minkowski spacetime. For weak data the whole spacetime evolving from 
the data has an asymptotic structure similar to Minkowski spacetime ||. Timelike infinity 
i + can be represented by a regular point on the grid. For strong data the spacetime will 
develop a singularity, there will not be any regular i + . In general, we expect strong field 
configurations to develop regions for which all future directed causal curves end in a singu- 
larity and the arising spacetimes will then have at least one event horizon. 
The next, slightly more complicated manifold is obtained by cutting out a ball from the 
interior of the ball as shown in figure 0. The interior £ has topology R x S 2 and there are 
two iS's evolving into two Js. This is the topology of a slice in the conformal picture of 
the Schwarzschild-Kruskal spacetime (see also figure ^). By Birkhoff's theorem spherically 
symmetric data with this topology of S will evolve into the Schwarzschild-Kruskal spacetime 
with its spacelike singularity in the future. For data close to Schwarzschild data there will 
probably be a spacelike singularity in the future and there will then be an event horizon 
around each J . Therefore it is justified to speak in the situation of figure of a one black 
hole spacetime, although there may be data for which no singularity appears. 
By cutting out N balls from the interior, spacetimes with N + 1 J~s are obtained. Figure |3| 
shows the N = 2 case. 
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FIG. 3. Data evolving into a space- 
time with an asymptotic structure as 
a two black hole spacetime with three 
disjoint null infinities. The third di- 
mension is again suppressed. The <S's 
are spheres. 



FIG. 4. Data evolving into a space- 
time with an asymptotic structure as 
the A3 spacetime. The third dimen- 
sion is again suppressed. The grid 
boundaries in the vertical dimension 
are identified as well as the boundaries 
in the suppressed dimension. The «S's 
are tori. 



Spacetimes usually interpreted as N black holes with N+l J\ have hyperboloidal hy- 
persurfaces like the one shown. Whether, in the case N = 2 for simplicity, an observer at J\ 
actually sees one or two black holes depends on the slice and the data. If the event horizon 
of J\ consists of two pieces, there are two separated black holes which may merge later, i. e. 
the two pieces of the event horizon merge. 



2. Ss with toroidal and other boundaries 



The examples discussed in the preceeding subsection are all asymptotically flat space- 
times which are only a subset of the asymptotically regular spacetimes. Except for one 
case spacetimes which are asymptotically regular, but not asymptotically flat, will not be 
discussed here. 

Although nowadays computers provide the resources for doing three dimensional calculations 
with reasonable resolution, high accuracy requirements as needed for certain interesting sce- 
narios may require calculations with higher resolution. Much higher resolutions can be 
achieved by analysing spacetimes with one Killing vector and they may also be easier to 
understand as they are easier to imagine and to visualise. In asymptotically flat spacetimes 
a Killing vector field with closed orbits vanishes at certain points and therefore there is an 
axis. Coordinate systems adapted to a one Killing vector symmetry become singular at the 
axis, the evolution equations or the variables are singular there, and numerical instabilities 
are very difficult to avoid. So it would be very interesting to have spacetimes available which 
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admit one Killing vector without having an axis, which possess null infinity, and which con- 
tain gravitational radiation. Those spacetimes are similar to the Schwarzschild solution, 
which is a spherically symmetric solution without a centre. 
This was the motivation of B. Schmidt to investigate these spacetimes 



15] . According to 



their simplest representative, the A3 solution in the classification of Ehlers and Kundt ||16|| , 
they are called asymptotically A3. Figure [| shows the initial slice for an asymptotically 
A3 spacetime. At the boundary of the suppressed z direction and the vertical y direction 
points are identified. Both <S's are tori. The A3 solution has a Killing vector with circular 
orbits but without axis. Many other asymptotically A3 solutions with two Killing vectors 
are known f|TJ,|l^]. In contrast to the spherically symmetric Schwarzschild solution many of 
them contain gravitational radiation ||18|| . 

Asymptotically A3 solutions provide an excellent testbed for numerical codes as they allow 
us the testing of radiation extraction procedures against exact solutions. They have been 

ITIJPI. Details about the tests of the 



intensely used for testing in the conformal codes 
author's code on these examples are going to be reported in [14 



C. J-fixing 

In the hyperboloidal initial value problem the boundaries S of S to are parts of ingoing null 
hypersurfaces. For a vanishing shift vector the coordinate area of S t shrinks, the physical 
part of the hypersurface is represented by less and less gridpoints. This loss of resolution 
is a wide-spread objection against the numerical application of conformal techniques. But 
this objection is unjustified as the shift vector can always be chosen in such a way that the 
location of S t is fixed and no gridpoints are lost. We call such a choice " j7"-fixing" . 
In the spherically symmetric case the author has used such a choice of coordinates in accuracy 
tests leading to 0. There were two reasons why this choice was not pursued further: Firstly, 
the accuracy obtained was lower than the one for the case of an "infalling" J even for late 
times and, secondly, i + was for that particular choice of coordinates no longer on the grid, 
the time direction was "decompactified" [unpublished work]. 

J. Frauendiener has generalised the j7-fixing JXTJ] to the non-spherically symmetric case. In 
the variables used in this paper the shift vector on J has to be chosen as follows: 

N a 1 7 = -^-h ab n b . (is) 

It is easy to see that (p~5|) plugged into equation ( |13h|) yields d t Vt\j = 0. The expression is 
regular on J as V a f2 is a non- vanishing null vector on J and therefore Qq 7^ 0, but Qq m ay 
become singular on M/S. As an alternative everywhere regular on M one can use: 

Since N a depends on variables, equation fll3c|) contains second derivatives of variables mak- 
ing the system of equations ( |i~3|) formally "parabolic" . But those derivatives can be elim- 
inated by the use of the other equations yielding a system first order in space and time 
again. J. Frauendiener has shown that the conformal field equations in spinor formulation 
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11]. He also demonstrated that ^-fixing 



stay symmetric hyperbolic even with ^-fixing 
works in his two dimensional code. 

By using a j7-fixing choice of shift there is no loss in freedom of specifying the shift vector 
on M, as with a j7-fixing shift iV a the new shift N a + QN a is also j7-fixing for arbritrary 
N a . 



IV. A WAY TO OBTAIN STABLE BOUNDARIES 

When solving symmetric hyperbolic equations numerically, the discretisation scheme 
used must be altered at the boundaries of the grid. It is not only necessary that the treat- 
ment must be consistent with the initial boundary value problem of the equations, the 
interior scheme and the boundary treatment must also fit. Otherwise gridmodes not related 
to a solution of the continuous equation arise from the boundaries and trigger instabilities. 
Even for very simple linear equations like the advection equation many discretisations of the 
boundary are unstable, although they look reasonable. 

In general relativity the problem is worse, even the analytic boundary value problem is not 
completely understood ||19|| . On a timelike boundary the evolution equations and the con- 



straints must be solved simultaneously. This consistency requirement determines the number 
of functions which can be specified freely. If the constraints are violated at the boundary, 
the numbers calculated at the gridpoints influenced by the boundary will not approximate 
any solution of the Einstein equation. The grid is filled with "invalid points" very quickly, 
independently of the grid spacing. To make things even worse the error made cannot be 
estimated by the size of the violation of the constraints HI. Therefore a discretisation must 
not only be stable but also consistent with constraint propagation. 

Even if one could solve this boundary problem there would still remain the problem of how 
to suppress artificial backscattering of an outgoing wave when it hits the boundary. Already 
for the wave equation in two and in three space dimensions that problem is not solvable by a 
local procedure [^T|, as the uniquely determined condition for "no backscattering" involves a 



non-local pseudo-differential operator. Stable local approximations to this non-local pseudo- 
differential operator show typically an artificial reflection of the order of percent | f2Tf . 
If it were necessary to calculate a solution of the conformal field equations on the whole grid, 
the situation would be similarly hopeless. But we do not have to calculate a solution of the 
conformal field equations on the whole grid. As already pointed out, MflP(S) = D(MflS). 
Therefore one can modify the equations and/or the data outside M U J . 
The "boundary treatment" we are going to present utilises this property of the conformal 
field equations and is independent of the exact form of the time evolution equations. To 
simplify the writing we continue by writing our time evolution equations in the general form 

dJ + A%l-b = 0. (17) 

The A* are quadratic matrices depending on /, b is a vector also depending on /. 

To obtain equations whose boundary treatment has already been analysed the equations ( |17|) 

are modified to 

dj + (1 - a(Q)) r%f_ + a(Q) (A%f ~ k) = 0, (18) 
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with a sufficiently often differentiable function which is for Q < Qq < ^1 < 

and 1 for O > Oi. r % is an outward pointing vector field. For Q > max grid boundary f2 one 
obtains near the boundary an equation with well-analysed boundary discretisations for many 
interior schemes. For an overview and the proofs of stability see [p2|p3| . The local order of 
convergence at the boundary may be one order lower than the convergence order n of the 



interior scheme without endangering the overall order n [|24 . 

The equations near the boundary are outward directed advection equations. Numerical noise 
created by a small transition zone [Qq, flj] is transported to the boundary. Since Q does 
not necessarily stay less than Q , the stability property of the boundary may get lost. By 
choosing r % = this can be avoided, the time evolution near the boundary gets frozen, one 
obtains: 

d t f + a{n)(A i d i f-b)=Q. (19) 

Of course the possibilities are not restricted to the cases discussed. One could e. g. try to 
combine the advantages of both methods described by only setting r l = in the equation 
for Q. 



The change of the system ( |TTD to (|IBD or ( |T9| ) will not change the solution on M. Nevertheless 
the numerically calculated numbers may change on M. In a numerical calculation we have 
to fulfil the Courant-Friedrich-Levy condition, hence the numerical range of influence of a 
point is larger than or equal to the analytic range of influence, i. e. values at gridpoints in 
M may numerically depend on values outside J . For a convergent scheme the dependence 
of the calculated solution on the values in the difference A of the two ranges must go to zero 
with the convergence rate. The error made by changing the data or the equations in A is 
of the same size, order, and nature as the discretisation error. This error is automatically 
taken into account when determining the error bars of the calculated values by a convergence 
analysis. Therefore, although the numbers actually calculated as an approximation to the 
solution may be changed by the boundary treatment, the physical predictions made from 
the numerical experiments do not change, since the physical predictions of a numerical 
experiment are the numbers modulo the error bars. 

Beside the simplification of the boundary treatment the change of the equations outside 
M has other positive effects. The equation ( |T7| ) may develop singularities everywhere in 
the conformal spacetime |3||§]. Since the equations of system (p~8|) or ( |19|) cannot develop 
singularities for Q < Q the appearance of a singularity outside M becomes less likely. 
The treatment of the grid boundaries described above is also independent of the choice 
of the gauge source functions and the coordinate representation of the light cone which 
depends on the gauge source functions. For eigenfeld methods, also applicable as boundary 
treatment [|TT|], this is not the case. The choice of the gauge source functions at the boundary 
influences the number of outward and inward propagating fields and, therefore, the boundary 
treatment. 



V. ASPECTS OF OUR GEOMETRICAL SETTING 

In the preceeding sections we have deduced a set of equations, which is equivalent to the 
Einstein equation and whose form is adapted to the requirements of numerics. To conclude 
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this paper, we will now describe on a heuristic level, based on the rigorous statements of the 
preceeding sections, how we proceed to calculate a spacetime from the initial slice. We do 
this on two examples of general interest and compare on those with alternative approaches. 

In the first example we illustrate by means of a Penrose diagram what part of a spacetime 
we can in principle calculate. To use a sufficiently rich example, i. e. an example for which 
inner and outer boundaries were believed to be needed, but still a simple enough to possess 
a depictable Penrose diagram, we have chosen the Schwarzschild-Kruskal spacetime. 
In the second example we explain on the scenario of an energy/matter concentration collaps- 
ing to a black hole how gravitational radiation propagates through the conformal spacetime. 

A. Evolving spacetime from an initial slice 

In figure [5] the Penrose diagram of the Schwarzschild-Kruskal spacetime is shown. 



singularity 




singularity 



FIG. 5. Penrose diagram of the Schwarzschild-Kruskal spacetime showing 
the region which can be calculated in the conformal approach. All points with 
the exception of the two spacelike infinities and (i°)2 represent an S 2 . 

In the numerical construction of this spacetime from an extended hyperboloidal initial 
value problem we give appropriate data on an extended hyperboloidal slice S to . The grid 
boundaries are placed outside the region representing the physical part of the initial slice. 
By solving the conformal field equations we calculate a slicing E$ of the future of Et . 
In principle there is no reason, assuming appropriate coordinates, which prevents us from 
calculating the whole domain of dependence of the physical part E to of E to up to the sin- 
gularity, i. e. to calculate the solution of the Einstein equation in the shaded region of the 
figure. The conformal field equations can be viewed as the Einstein equation for an artificial 
spacetime with an artificial matter field, namely the rescaling factor. With this analogy in 
mind it becomes obvious that the relation between coordinate choice and exhaustive slicing 
of the spacetime is similar in conformal and physical spacetime. 

We now describe how the same situation is approached in two alternative numerical 
approaches, namely the approach of solving the initial boundary value problem in physical 
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spacetime and the so-called Cauchy characteristic matching approach. 

Figure |6] shows the Penrose diagram of the first case, the approach by solving the initial 

boundary value problem in physical spacetime as described in f25| . 



singularity 




FIG. 6. The part of the Penrose diagram 
of Schwarzschild spacetime which is calculated 
when solving the initial boundary value prob- 
lem in physical spacetime. Again every point 
represents an S 2 . 



Data are given on a compact part of a spacelike initial slice, which is covered by grid- 
points, and on the boundaries. To be representative for the whole slice, which has topology 
R x S 2 , the compact part should have topology I x S 2 , where / is a closed interval. The 
inner edge of that closed interval / is placed somewhere inside the event horizon, since then 
its treatment will not influence the physical predictions. For technical reasons the inner 
"boundary" is often placed near the apparent horizon and therefore called apparent horizon 
boundary, although the apparent horizon is a null or spacelike hypersurface. In contrast to 
the inner edge the outer edge of the grid is treated in such a way that it forms a timelike 
boundary under time evolution. Due to the problem of not having a consistent outer bound- 
ary treatment without artificial back scattering, the numerical solution will be consistent 
with the Einstein equation in the hatched triangle only. The actual size of the hatched 
region depends on how much gridpoints one can afford to provide to cover the parts of the 
hypersurface far away. Assuming a consistent outer boundary treatment one would get the 
shaded region. Null infinity can never be covered, even after rescaling, for the spacelike 
surfaces used go to spacelike infinity as sketched by the dashed lines. 

The approximate radiation extraction methods suggested require being "far away" to give 
a reasonable approximation for the actually emitted radiation. Nevertheless reading off ra- 
diation at a finite distance causes an approximation error. Hence the error made will not 
converge to zero when refining the grid. It has therefore been suggested to combine the grid 
refinements with placing the outer boundary further out. There are two problems with this 
"solution" . Firstly, the number of required gridpoints grows much faster than in the case of 
only refining and this additional growth certainly outweighs the benefits from computational 
science efforts to optimise memory footprint and execution speed per gridpoint. Secondly, 
whenever we place the boundary at a different distance we solve a different problem. This 
change of the problem does not occur in the conformal method, since the grid contains for 
every refinement level the conformally completed spacetime and null infinity. 
The second alternative we compare with, Cauchy characteristic matching ||26|| , is shown in 
figure 0. 
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FIG. 7. The part of the Penrose diagram of 
Schwarzschild spacetimes which is calculated in 
the Cauchy characteristic approach. Again ev- 
ery point represents an S 2 . 



Near the outer edge and near the apparent horizons the spacetime is sliced by null hyper- 
surfaces. The inner null hypersurfaces are cut off inside the apparent horizon. The region 
between the null hypersurfaces is evolved via a Cauchy problem. As null hypersurfaces un- 
avoidably develop caustics for strong data, the Cauchy slicing in the middle is necessary and 
must stretch out sufficiently far. Matching the two schemes introduces artificial boundaries 
and causes stability problems f26|| . Although null infinity is mapped to gridpoints by a 
coordinate rescaling, null infinity remains a true numerical boundary as the equations used 
are singular there. 

As a conclusion we find that in all approaches only part of the whole spacetime can be 
calculated. But this does not favour the other approaches discussed, because none of them 
allows us to calculate a larger part of the Penrose diagram. Even when using the data on 
S to to integrate backwards in time, as done in [figure Va], the neighbourhood of the i°s 
cannot be obtained. To obtain the whole Penrose diagram one would have to start with 
an initial slice containing spacelike infinity. The conformal metric at i° is in general not 
sufficiently smooth, some of the variables used in the conformal field equations are not well- 
behaved at i°. Therefore the more advanced and complicated conformal techniques of |2£ 
would be required. 



B. "Tracking radiation sufficiently far" or "The need for apparent horizons" 

Probably the simplest scenario which involves gravitational radiation and for which ap- 
parent horizon boundary conditions have been regarded as necessary is an energy/matter 
concentration collapsing to a black hole. 

Figure |8] shows in a diagram how this scenario is calculated in the conformal approach. Its 
correctness has been justified by the model calculations discussed in 
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FIG. 8. Cauchy slices in conformal FIG. 9. Hyperboloidal slices in physical 

spacetime (= hyperboloidal slices in spacetime 
physical spacetime) 



In the situation drawn a singularity avoiding slicing of the conformal spacetime is in no 
way pathologic — null infinity, timelike infinity, and the singularity are all at finite confor- 
mal coordinate time. A pulse of radiation will propagate outward with a slope of order one 
and will cross J at finite conformal coordinate time. When this pulse crosses J we can read 
off the radiation which is seen by an observer at infinite physical distance. Although the 
physical distance between the gridpoints becomes very large near J , the number of grid- 
points used for resolving a signal, as indicated by the width of the pulse, stays approximately 
constant. The mathematical reason behind this feature is the invariance of the light cone 
under conformal rescalings. 

If we translate figure [8] to physical spacetime coordinatized by physical time t and by physi- 
cal distance r we get figure H The hyperboloidal character of the hypersurfaces £f becomes 
obvious, for large distances the hypersurfaces approach null hypersurfaces. Since the hy- 
persurfaces are almost null far out, the finite, but far distance at which we want to read 
off radiation does not significantly influence the integration time required. In the example 
shown we would need to calculate a numerical time evolution for t = 35 wadm to be able to 
read off radiation at f = 100 ^adm- 

We also see that these "physical" coordinates do not very well represent the causal structure 
of the spacetime near and inside the event horizon — the singularity and the event horizon 
look like timelike lines. Furthermore we notice that in the interior the slices cease to be 
regular after a finite time, whereas in the exterior they exist forever. Therefore, if we would 
like to calculate the longtime behaviour as seen by an observer outside the event horizon we 
must either distort the hypersurface somewhere or cut out the interior part. 
Figure |IT] shows the situation in the physical spacetime with singularity avoiding slices, the 
slicing which is believed to be the most appropriate way for approaching the problem in 
physical spacetime. 
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FIG. 10. Cauchy slices in physi- 
cal space-time 



It is obvious that if we want to read off radiation at f = 100 h-adm? we need to calculate 
the time evolution up to t = 100 w^adm? an d if we want to read it off at f = 1000mADM 5 
we even need to integrate to t — 1000 mADM- Therefore, even if we are only interested in 
gravitational radiation but not in the longtime behaviour of the spacetime, it is necessary 
to numerically integrate for a long time, at least in the exterior. To not hit the singularity 
when doing so one again either has to deal with a significant distortion of the hypersurfaces, 
i. e. large gradients of the variables somewhere on the hypersurfaces, or one has to cut out 
the interior of the hypersurfaces, which is the purpose of using apparent horizon boundary 
conditions. 

From what has been said about figure |9| and |10| one might be tempted to utilize hyperboloidal 
slices in the physical spacetime for the numerical calculation. This would, unfortunately, 
introduce new problems. They would arise from the fact that the evolution equations require 
spacelike surfaces but the hyperboloidal slices are almost null in the far zone. Technically, the 
lapse and/or shift degenerate when a hypersurface becomes almost null. Since in conformal 
spacetime the hyperboloidal hypersurface is a normal Cauchy surface even at null infinity 
and beyond, this problem does not arise in the conformal approach. 



VI. CONCLUSION 

Starting from the conformal field equations we have derived a symmetric hyperbolic first 
order system determining the time evolution. The system obtained is in a form which allows 
us to apply standard techniques to discretise in the interior as well as on the grid boundaries. 
The stability of the resulting schemes is then guaranteed by theorems in the literature. 
We have seen that an application of only these techniques allows us in principle to calculate 
the complete future of the initial data for many interesting scenarios, e. g. N black holes. 
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